4. Modulazione di ampiezza

signal
filter
spectrum
Autore/Autrice

Mariolino De Cecco

Data di Pubblicazione

24 febbraio 2025

Sommario

Esercizi sulla modulazione di ampiezza delle misure mediante analisi dello spettro e filtraggio offline.

1 Librerie

Quanto segue richiede questi pacchetti:

library(tidyverse)
library(plotly)
library(gsignal)

Se necessario installarli con install.packages() o mediante la GUI.

2 Funzioni di utilità

Per comodità definiamo alcune funzioni di utilità.

Generazione di un segnale sinusoidale con diverse armoniche, come definite in una tabella (pars) con le colonne w (ampiezza), f (frequenza) e phi (fase) per le varie frequenze (una riga per ogni armonica):

signal <- function(t, pars, rad = FALSE) { 
  stopifnot(is.data.frame(pars))
  with(pars, {
    if (!rad) {
      phi <- phi/180*pi
      f <- 2*pi*f
    }
    map_dbl(t, \(t) sum( map_vec(seq_along(w) , \(i) w[i]*sin(t*f[i] + phi[i] ))))
  })
}

3 Simulazione del segnale

Costruiamo l’uscita di uno strumento di misura, la sua modulazione in ampiezza, il rumore interferente e grafichiamo gli spettri.

Una nota su plotly

Nel corso di statistica viene spiegato l’uso di GGPlot2. Quest’ultimo è studiato per l’efficacia e chiarezza di presentazione, per le quali è superiore a plotly.

È comunque possibile rendere i grafici GGPlot2 interattivi, passandoli alla funzione ggplotly().

Si noti tuttavia che avere molti grafici interattivi (e con molti punti) rallenta parecchio il caricamento della pagina nel browser: meglio evitarlo per i grafici “intermedi” e riservarlo ai grafici finali.

Definiamo alcuni parametri generici e la tabella delle armoniche:

# Parametri
fc <- 2000        # Frequenza di campionamento (Hz)
fm <- 10          # Frequenza del segnale di misura (Hz)
fr <- 15          # Frequenza del rumore (Hz)
fp <- 300         # Frequenza della portante (Hz)
m <- 0.8          # Indice di modulazione (0 < m < 1)
duration <- 10    # Durata del segnale (secondi)

# tabella delle armoniche del segnale di misura
pars_m <- tibble(
  w = c(1, 0.5, 0.3), 
  f = c(fm, 2*fm, 4*fm), 
  phi = c(0, 0, 0)
)

pars_m %>% knitr::kable()
w f phi
1.0 10 0
0.5 20 0
0.3 40 0

Generiamo il segnale di misura, perturbato da un rumore normale:

ym <- tibble(
  t = seq(0, duration, by = 1/fc),
  y = signal(t, pars_m),
  yn = y + rnorm(length(t), 0, pars_m$w[1]/10)
) %>% 
  mutate(
    f = 0:(length(t)-1)/max(t),
    fft = fft(y)
  ) 

ym %>%  
  mutate(intensity = Mod(fft) / length(t)*2,) %>% 
  ggplot(aes(x=f, y=intensity)) +
  geom_line()

Generiamo il rumore interferente con il suo spettro. Parametri delle armoniche corrispondenti al rumore interferente:

pars_r <- tibble(
  w = c(1, 0.5, 0.8), 
  f = c(fr, 3*fr, 9*fr), 
  phi = c(0, 40, 90)
)
pars_r %>% knitr::kable()
w f phi
1.0 15 0
0.5 45 40
0.8 135 90

Aggiungo il rumore interferente alla tabella ym:

ym <- ym %>% 
  mutate(
    y_in = signal(t, pars_r),
    fft_in = fft(y_in)
  )

Ottengo il grafico dello spettro:

ym %>% 
  select(f, fft, fft_in) %>% 
  pivot_longer(
    contains("fft"), 
    names_to = "signal",
    values_to = "fft",
    names_transform = ~ ifelse(.=="fft", "signal", "intf. noise")
  ) %>% 
  mutate(
    intensity = Mod(fft) / length(t)*2
  ) %>% 
  ggplot(aes(x=f, y=intensity, color=signal)) +
  geom_line()

Tidy data

Cerchiamo più possibile di utilizzare dati tidy, cioè una variabile per colonna, un’osservazione per riga. In questo modo realizzare grafici con serie multiple è particolarmente rapido (basta identificare la variabile di classificazione, cioè il nome della colonna che contiene la chiave di raggruppamento).

Per questo motivio, qui e di seguito useremo pivot_longer.

Aggiungo anche la modulazione:

# Aggiungo una mappa da nomi delle colonne -> nomi in legenda
fft_names=c(
  "fft" = "signal", 
  "fft_in" = "intf. noise", 
  "fft_m" = "modulation"
)

# Aggiungo la modulazione
ym <- ym %>% 
  mutate(
    y_m = (1 + m * yn) * cos(2 * pi * fp * t),
    fft_m = fft(y_m)
  )
# Grafico dello spettro
(ym %>% 
  select(f, contains("fft")) %>% 
  pivot_longer(
    -f, 
    names_to = "signal",
    values_to = "fft",
    names_transform = ~ fft_names[.]
  ) %>%
  mutate(
    intensity = Mod(fft) / length(t)*2
  ) %>% 
  ggplot(aes(x=f, y=intensity, color=signal)) +
  geom_line() +
  xlim(c(0, 500))) %>% 
  ggplotly()
Alternativa SENZA pivot_longer

Ridotto all’osso, l’ultimo grafico è (evito ggplotly per brevità):

ym %>% 
  select(f, contains("fft")) %>% 
  pivot_longer(-f) %>%
  mutate(
    intensity = Mod(value) / length(t)*2
  ) %>% 
  ggplot(aes(x=f, y=intensity, color=name)) +
  geom_line()

Se non voglio usare pivot_longer() devo aggiungere tre colonne, una alla volta (stando attento a ripetere la formula correttamente: se scrivessi fft_m = Mod(fft_n) / length(t) * 2 otterrei un errore molto difficile da individuare!), con il modulo della FFT e nel plot devo aggiungere un layer alla volta (stando asttento a specificare correttamente il nome della serie come variabile color):

ym %>% 
  select(f, contains("fft")) %>% 
  mutate(
    fft = Mod(fft) / length(t) * 2,
    fft_in = Mod(fft_in) / length(t) * 2,
    fft_m = Mod(fft_m) / length(t) * 2
  ) %>% 
  ggplot(aes(x=f)) + 
  geom_line(aes(y=fft, color="fft")) +
  geom_line(aes(y=fft_in, color="fft_in")) +
  geom_line(aes(y=fft_m, color="fft_m"))

4 Effetto interferente sul segnale modulato

Adesso che abbiamo il segnale modulato separato in frequenza dal rumore possiamo ‘esporlo’ al rumore interferente.

Osserviamo prima la modulazione:

ym %>% 
  ggplot(aes(x=t, y=y_m)) + 
  geom_line()

Combiniamo segnale e modulazione:

ym <- ym %>% 
  mutate(
    y_m_in = y_m + y_in,
    fft_m_in = fft(y_m_in)
  )

Ora ym è così composta:

ym %>% head() %>% knitr::kable()
t y yn f fft y_in fft_in y_m fft_m y_m_in fft_m_in
0.0000 0.0000000 0.1418529 0.0 0.00e+00+0.0000000i 1.1213938 1.121394+0.0000000i 1.1134823 3.829540+0.000000i 2.2348761 4.950934+0.000000i
0.0005 0.1004060 0.1040175 0.1 2.10e-06-0.0132478i 1.1483848 1.121397-0.0073390i 0.6366972 4.564120+5.526609i 1.7850820 5.685517+5.519270i
0.0010 0.2000641 -0.0226319 0.2 8.30e-06-0.0265019i 1.0386500 1.121407-0.0146797i -0.3034221 -3.960835-5.672208i 0.7352280 -2.839429-5.686887i
0.0015 0.2982363 0.3224292 0.3 1.87e-05-0.0397689i 0.8266720 1.121422-0.0220241i -1.1963752 4.411363-12.118149i -0.3697032 5.532785-12.140174i
0.0020 0.3942043 0.4672451 0.4 3.33e-05-0.0530552i 0.5637101 1.121445-0.0293739i -1.1114244 -5.362537+4.257905i -0.5477143 -4.241092+4.228531i
0.0025 0.4872785 0.5019278 0.5 5.21e-05-0.0663672i 0.3085893 1.121474-0.0367308i 0.0000000 4.785527+5.336739i 0.3085893 5.907001+5.300008i
ym %>% 
  select(t, y_m, y_m_in) %>% 
  pivot_longer(-t, names_to = "series", values_to = "value") %>% 
  ggplot(aes(x=t, y=value, color=series)) +
  geom_line() +
  facet_wrap(~series, nrow=2) +
  xlim(c(0, 0.5))

Osserviamo lo spettro:

fft_names["fft_m_in"] <- "mod. + intf."
ym %>%
  select(f, contains("fft")) %>% 
  pivot_longer(
    -f, 
    names_to = "signal",
    values_to = "fft",
    names_transform = ~ fft_names[.]
  ) %>% 
  mutate(
    intensity = Mod(fft) / length(t)*2
  ) %>% 
  ggplot(aes(x=f, y=intensity, color=signal)) +
  geom_line() +
  xlim(c(0,500))

Esercizio

Estrarre il segnale in uscita originale a partire da quello affetto da rumore bianco, modulato ed affetto da rumore interfente ovvero da ym$y_m_in

4.1 In un colpo solo

Nella parte precedente, la tabella ym è stata costruita gradualmente, allo scopo di illustrare la composizione del segnale simulato. Guardando le operazioni nel loro complesso, la procedura è in realtà più semplice:

tibble(
  t = seq(0, duration, by = 1/fc),
  y = signal(t, pars_m) + 
      rnorm(length(t), 0, pars_m$w[1]/10),    # segnale 
  fft = fft(y),                               # e sua FFT
  y_in = signal(t, pars_r),                   # rumore interferente
  fft_in = fft(y_in),                         # e sua FFT
  y_m = (1 + m * y) * cos(2 * pi * fp * t),   # modulante
  fft_m = fft(y_m),                           # e sua FFT
  y_m_in = y_m + y_in,                        # modulante + interferente
  fft_m_in = fft(y_m_in)                      # e sua FFT
) %>%
  # aggiungo la frequenza
  mutate(
    f = 0:(length(t)-1)/max(t),
  ) %>% 
  # prendo solo le colonne frequenza e fft
  select(f, contains("fft")) %>% 
  # riorganizzo la tabella
  pivot_longer(
    -f, 
    names_to = "signal",
    values_to = "fft",
    names_transform = ~ fft_names[.]
  ) %>% 
  # aggiungo l'intensità
  mutate(
    intensity = Mod(fft) / length(t)*2
  ) %>% 
  ggplot(aes(x=f, y=intensity, color=signal)) +
  geom_line() +
  xlim(c(0,500))

5 OMETTERE E FAR FARE COME ESERCIZIO

Passo 1: Filtraggio e Demodulazione

Estraiamo il segnale modulato mediante un filtro passa-alto. Di conseguenza è rimasto solo il segnale modulato che va quindi demodulato per ottenere una stima del segnale in uscita originario

cutoff <- 2 * max(pars_r$f) / fc  # Frequenza di taglio (2 volte la massima frequenza del rumore)
bf <- butter(4, cutoff, type = "high")
ym <- ym %>% 
  mutate(
    y_m_f = filtfilt(bf, y_m_in) # filtfilt non introduce ritardo
  )

Passo 2: Demodulazione

Moltiplicazione per la portante:

ym <- ym %>% 
  mutate(
    demod_raw = 2 * y_m_f * cos(2 * pi * fp * t) - 1
  )

Passo 3: Filtraggio Passa-Basso con filtro Butterworth

cutoff <- 2 * max(pars_m$f) / fc  # Frequenza di taglio (2 volte la massima frequenza del segnale informativo)
bf <- butter(4, cutoff, type = "low")
ym <- ym %>% 
  mutate(
    demod = filtfilt(bf, demod_raw)
  )

Confronto dei segnali:

(ym %>% 
  mutate(y_noisy = y + y_in) %>% 
  select(t, y, y_noisy, demod) %>% 
  pivot_longer(-t) %>% 
  ggplot(aes(x=t, y=value, color=name)) + 
  geom_line() +
  coord_cartesian(xlim=c(0, 0.5))) %>% 
  ggplotly()

Si noti come il segnale originario fosse affetto da rumore bianco, poi modulato e quindi sovrapposto ad un effetto interferente a tre armoniche. Alla fine del processo abbiamo ottenuto il segnale molto simile al segnale originario con una riduzione vistosa anche del rumore interferente